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We show that in the collision of two superfluid fermionic atomic clouds one observes the formation 
of quantum shock waves as discontinuities in the number density and collective flow velocity. Domain 
walls, which are topological excitations of the superfluid order parameter, are also generated and 
exhibit abrupt phase changes by tt and slower motion than the shock waves. The domain walls are 
distinct from the gray soliton train or number density ripples formed in the wake of the shock waves 
and observed in the collisions of superfluid bosonic atomic clouds. Domain walls with opposite phase 
jumps appear to collide elastically. 
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Quantum shock waves (QSWs) and solitons have been 
observed in a dilute atomic Bose gas a decade ago pQ, a 
result which generated a flurry of experimental and the- 
oretical work [2]-[5]. In classical hydrodynamics, shock 
waves appear as discontinuities in density, flow velocity, 
temperature and other characteristics of the fluid. In 
a fluid in motion one can observe another remarkable 
structure, solitons. Solitons are usually a manifestation 
of the competition between the dispersive mechanisms 
and nonlinearities at work in a fluid. While dispersion 
effects tend to smooth out discontinuities, nonlinearities 
sometimes oppose this effect instead, stabilizing a large 
local density variation [6 that leads to the appearance 
of a bright, dark or a gray soliton. The quantum fluid 
dynamics of a superfluid dilute atomic Bose gas is typi- 
cally described using the nonlinear Gross-Pitaevskii (gp) 
equation, and at sufficiently low temperatures there does 
not appear to be any need for dissipative processes. The 
accurate treatment of the real-time dynamics of a super- 
fluid Unitary Fermi Gas (ufg) requires a more complex 
approach, using an extension of the Density Functional 
Theory (dft) to superfluid systems and time-dependent 
phenomena, see Refs. [7HTT]. 

Recently, a new experiment reported on the observa- 
tion of QSWs in a UFG [12], see Fig. [I] The UFG is a 
system exactly in the middle of the BCS-BEC crossover 
and as such its properties often qualitatively interpolate 
between those of a Fermi and a Bose superfluid. However, 
while BCS and BEC systems are weakly interacting, the 
UFG is a strongly interacting superfluid in which the crit- 
ical temperature and the Landau critical velocity attain 
their highest values (in appropriate units). The veloc- 
ity equation used in Ref. [12 to model shock waves in 
the collisions of two UFG clouds neglects the "quantum 
pressure term" and a phenomenological viscosity term is 
added: 



in which one can choose the mass of the atom m = 1, v k 
and Vfc are the cartesian coordinates of v and V respec- 
tively, and [V0^| = \7 k vi + Viv k -2/35 k iV -v. Above 
ji[n] is the chemical potential in homogeneous matter at 
a given number density n, and we have suppressed the 
explicit dependence of the number density n(r,t), veloc- 
ity v(r,t) and external trapping potential V ex t(r : t) on 
space-time coordinates. In the superfluid phase both the 
shear viscosity and the "quantum pressure" terms are 
proportional to h and a priori it is not obvious that one 
can neglect one term, but not the other. Unlike the case 
of dilute Bose gases, where the QSWs were interpreted 
as "pure" dispersive shock waves [3] with no need for 
dissipative effects, the results of the experiment [12] on 
UFG received an interpretation similar to classical shock 
waves, in which dissipation plays a crucial role in the for- 
mation of the shock wave front. These two distinct inter- 
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FIG. 1. (Color online) The space-time evolution of the ID 
density profiles along the trap axis as reported in Fig. 2 of 
Ref. [12]. We have outlined part of the shock wave front with 
a white oval here and in Figs. [2] [4] [5] and [6] 
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pretations of the experiments on Bose and Fermi dilute 
gases are difficult to reconcile. In the BEC regime, the role 
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of dissipation is negligible (at least at the phenomenolog- 
ical level) and the QSW and the density ripples identified 
with soliton trains can be described by dispersive effects 
alone. Viscosity was introduced phenomenologically in 
Ref. [12 j, to a large extent in order to avert the onset of a 
"gradient catastrophe" [15] . At the same time, the UFG is 
widely accepted as a prime example of an almost perfect 
fluid [14 . In the UFG the bulk viscosity vanishes and the 
shear viscosity is at a minimum as a function of the cou- 
pling constant across the BEC-BCS crossover, see Refs. 
[T4HT6] and references therein, and thus we see no com- 
pelling theoretical arguments to include it in the present 
study. In this respect, the present approach is similar to 
other studies of dilute Fermi gases [TTHST] . Another sig- 
nificant limitation of the hydrodynamic approach, which 
was used in Ref. [12] , is its inability to describe quantum 
topological excitations (quantized vortices, domain walls 
(dw)), both of which have been observed in the similar 
experiments with bosons [IHS]. 

In the case of colliding UFG clouds, we observe the gen- 
eration of both QSWs and DWs. The DW excitation has 
been suggested in other simulations [T7H20] . The DWs are 
excitations of the superfluid order parameter and not the 
number density ripples identified as soliton trains trailing 
the wake of the shock waves, as discussed in Refs. [THS]. 
We will make this distinction in order to avoid confusion. 
We show that the number density of two colliding UFG 
clouds shows a behavior very similar to the one observed 
in experiment [T2] , In the wake of the QSWs we observe 
in addition the formation of DWs. The DWs emerge as 
quite sharp changes in the phase of the superfluid order 
parameter by 7r, and are correlated with minima of the 
number density. DWs propagate through the system at 
slower speeds than the QSWs and are topological excita- 
tions similar to quantum vortices. DWs always appear in 
pairs with opposite jumps of the order parameter phase 
and appear to collide essentially elastically with one an- 
other and with the system boundary. These phenomena 
are observed in the absence of any dissipation, which is 
expected to play a negligible role at temperatures close 
to absolute zero, see Ref. [I5j [16] and references therein. 

The extension of the dft to superfluid fermionic sys- 
tems and time-dependent phenomena has been described 
and applied to a number of phenomena in nuclear physics 
and the physics of cold gases [9HTT]. This approach 
is known as the Superfluid Local Density Approxima- 
tion (slda) and tdslda for its time-dependent ver- 
sion respectively, based on the simplest possible (un- 
renormalized) energy density functional: 
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and 



v c = <E v^u n are the kinetic, number and anoma- 
lous densities respectively, and V ex t is an external one- 
body potential. The dynamical evolution of the system is 



described by Bogoliubov-de Gennes-like equations (BdG) 
for the quasiparticle wave functions (qpwfs) {u ni v n ) 
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where the single-particle Hamiltonian, /i, and pairing po- 
tentials, A, are obtained by taking the appropriate func- 
tional derivatives of the energy density functional £ . The 
dimensionless constants a, f3 and 7 are fixed by the en- 
ergy per particle, pairing gap and quasiparticle spectrum 
obtained from Quantum Monte Carlo (qmc) calculations 
of the homogeneous infinite system. For the descrip- 
tion of the renormalization procedure, various technical 
details, and numerical implementation details see Refs. 
[TO] [TT] . Within SLDA various properties of a UFG are re- 
produced within a few percent accuracy, mainly limited 
so far by the accuracy of the QMC values for the pairing 
gap and for the effective mass [10] . Dimensional and sym- 
metry arguments, renormalizability, and Galilean invari- 
ance define tdslda uniquely and with an overall accu- 
racy not worse than 10%. Unlike the BdG approximation, 
in which interaction effects vanish in the absence of pair- 
ing correlations, the energy of both superfluid and normal 
phases at zero temperature are described accurately. The 
UFG has a condensation energy « 20% of the total inter- 
action energy. In the BdG approximation, however, the 
condensation energy is equal to the interaction energy in 
the case of an UFG. Even at unitarity, the meanfield en- 
ergy dominates the pairing energy, and that affects the 
dynamics accordingly. An UFG initially in a superfluid 
phase that is subjected to an external time-dependent 
agent can become normal and the Cooper pairs could be 
destroyed, but the particles will still be strongly inter- 
acting. In contradistinction, in the BdG approximation, 
which was used in studies of solitons [T7H20] , the normal 
phase is simulated as a non- interacting Fermi gas [TT] , 

We have performed simulations of the cold atom 
cloud collisions assuming that qpwfs have the struc- 
ture u n {r,t) exp(ik nz z)u n (x,y,t),v n (r,t) => 
ex.-p(ik nz z)v n (x,y,t) with periodic boundary conditions 
in the z-direction and a rather stiff harmonic confining 
potential in the ^-direction. The time-dependent trap- 
ping potential along the collision x-axis had a similar 
profile and time dependence as in experiment [12J. The 
solitons and the shock waves now are 2D in character, 
their stability properties are slightly different than in 3D, 
and the sound velocity is modified [22]. Typical results 
of these simulations are shown in Figs. [2j [3j [4j [5| and 
[6] One can zoom into the online figures in order to see 
details, which otherwise would escape the naked eye. 

The simulation results, Fig. [2] show remarkable 
similarities with the experiment, Fig. [l] The QSWs 
speed in experiment (v qsw /vF ~ 0.35) and simulations 
(vqsw/vf ~ 0.25) agree within w 25 - 30%. The dif- 
ferences can be ascribed to various experimental un- 
certainties (in particular the particle number), different 
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FIG. 2. (Color online) The space-time evolution of the num- 
ber density profile n(x,0,£) along the collision axis in TD- 
SLDA simulation of the collision of two UFG clouds. Here Jzf 
and sf = kp/2 (h — m — 1) are the initial values of the Fermi 
wave vector and energy at the center of clouds. The region 
where one elastic collision of two dws occurs is outlined with 
a white circle here and in Figs. [4] [5] and [5] 




FIG. 3. (Color online) Three consecutive frames showing the 
absolute magnitude of the pairing field \A(x,y,t)\ in the xy- 
plane at times teF — 30, 350, 690, see also Figs. [2]and[4] The 
x- and ^/-directions (not shown to scale here) have an aspect 
ratio of ~ 30. By the time of the second frame the QSWs 
have emerged from the cloud, dws appear as planar number 
density depletions with a width comparable to the diameter 
of a quantum vortex, and can be highly visible using standard 
techniques [23 . We have observed the formation of dws so 
far only in traps with elongations larger than in experiment 
|12j . In simulations of clouds with smaller aspect ratios we 
observe only QSWs with density profiles very similar to Ref. 

trap shapes (Vtrap(x,y, z) oc uo 2 x 2 + uo^y 2 + z 2 ) and 
Vtrap(%,y, z) oc u 2 x 2 + ujj_y 2 in experiment and simula- 
tions respectively) and sizes, different cloud aspect ratios, 
and also different amount of initial collisional kinetic en- 
ergy injected into the system. In spite of being confined 
in the ^-direction in a harmonic potential, the DW are 
planes perpendicular to the collision x-axis, see Fig. [3] 
However, in Ref. [12 the experimental set-up prevented 
the authors from observing the DWs. The images cor- 
responding to various frames reported there were taken 
in different realizations of the two colliding clouds. On 
one hand, the phase differences of the two initially sepa- 
rated condensates are likely random and maybe difficult 
to control, see Ref. [24]. On the other hand, the den- 




FIG. 4. (Color online) The magnitude of the pairing gap 
|A(x,0,£)| corresponding to Fig. [5] The dws appear as sig- 
nificant depletions of the pairing field, and always appear in 
pairs with opposite phase jumps, see also Figs. [5] and [6] 

sity profile fluctuations from shot-to-shot in Ref. [12] 
point to a rather low spatial resolution attained in these 
measurements (see Fig. 2 in Ref. [12 J, which explains 
why DWs have not been observed in this experiment. We 
have performed simulations by varying the initial rela- 
tive phase of the condensates. While the overall picture 
of the collisions remains unchanged, the number of DWs 
created varies. The density ripples in the wake of the 
shock waves discussed in experiments with Bose dilute 
clouds [THS] and interpreted there as a soliton train, are 
formed here as well. By zooming in the online Figs. [2] 
and [4] one can notice that before the shock wave is formed 
well defined matter wave interference occurs. The discon- 
tinuity in the number density and order parameter at the 
wake of the shock wave is accompanied by a similar dis- 
continuity in the collective flow velocity, see Fig. [6] The 
DWs which form in the wake of the shock wave have lower 
speeds. The width of the domain wall density depletion 
is comparable with the diameter of the density depletion 
of the quantized vortex core in a UFG [7J [8] . This as- 
pect might present a challenge if one were to attempt a 
direct observation. However, the experimental technique 
implemented to visualize the quantum vortex lattice in 
a UFG [23] can be implemented to put in evidence DWs, 
particularly since a spatial imprint left by a single DWs 
is a much larger 2D-structure, as compared to the ID 
imprint left by vortex core, see Fig. [3] In sufficiently 
wide traps DWs might develop "snake" instabilities as in 
Ref. [1 , which would be a clear "smoking gun" of DWs 
formation. 

In summary, we have shown that QSWs and DWs are 
formed in the collision of two initially independent super- 
fluid fermionic clouds. The QSWs manifest themselves as 
rather sharp discontinuities in the number density, super- 
fluid order parameter and collective flow velocity. They 
reflect essentially elastically from system boundaries, but 
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FIG. 5. (Color online) The phase of the pairing gap 
arg A(x, 0, t)/it. 
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FIG. 6. (Color online) The x-component of the collective 
flow velocity field v x (x, 0, t) along the axis of collision. At the 
front of the two shock waves the velocity field undergoes a 
rapid change, and the matter flows in opposite directions. 



with significant reduction in their intensity. The shock 
waves lead to the formation of DWs, topological excita- 
tions of the superfluid order parameter similar to quan- 
tum vortices. As in the case of vortices, DWs show-up 
as significant depletions of the number density, but over 
much wider spatial regions [TJ [23] . The phase of the 
superfluid order parameter changes quite abruptly by 7r 
across a domain wall, and DWs appear always in pairs 
propagating in opposite directions. The DWs reflect elas- 
tically from the system boundaries. However, they even- 
tually dissipate, in particular when colliding with one 
another, if they have the same phase jump. DWs with 
similar phase jumps can propagate at different speeds 
and when they catch one another often annihilate. DWs 
with opposite phase jumps appear to collide elastically. 
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